OpenEO Platform - User Story 3

On-demand processing of Sentinel-1 data:

• ARD compliant
• With custom parametrization

On-demand processing of Sentinel-1 data

Compare the data given the different options: ARD metadata, terrain correction?, sigma, gamma

• Visualize the metadata
• Sigma to Gamma Ratio
• Scaling from linear to decibel?

OpenEO conversion to dB in the process graph.

• Comparison in Bolzano area
• Show that changing a parameter the output changes.
• Show sigma to gamma ratio as another raster.

from eo_utils import *

Select the area of interest in the next map

center = [46.49, 11.35]
zoom = 12

eoMap = openeoMap(center,zoom)
eoMap.map

openEO graph creation and processing for S1 data

1. Open the connection with VITO back-end. Please use the openeo-auth tool to configure credentials. https://open-eo.github.io/openeo-python-client/auth.html#config-files-and-openeo-auth-helper-tool

connection = openeo.connect("https://openeo-dev.vito.be").authenticate_basic()

1. Define range of interest in time from Sentinel-1

if(len(eoMap.bbox) == 0):
    mapBox =eoMap.map.bounds     
    bbox = [ mapBox[0][1],mapBox[0][0],mapBox[1][1],mapBox[1][0]]
else:
    bbox = eoMap.bbox
bbox

[11.236267089843752, 46.46825122515853, 11.457710266113283, 46.51552383029837]

collection      = 'SENTINEL1_GRD'
spatial_extent  = {'west':bbox[0],'east':bbox[2],'south':bbox[1],'north':bbox[3],'crs':'EPSG:4326'}
temporal_extent = ["2021-01-01", "2021-01-10"]
bands           = ["VV","VH"]

s1 = connection.load_collection(collection,spatial_extent=spatial_extent,bands=bands,temporal_extent=temporal_extent)

1. Apply the ARD compliant SAR processing.

s1bs = s1.ard_normalized_radar_backscatter()

1. Visualize the openEO process graph.

s1bs.to_graphviz()

1. Create a new batch job on the back-end and ask to process it.

If our area of interest is small, we can also do a direct request, but this will not return the CARD4L json metadata.

#skip evaluating the rdd, this is a small performance improvement
#s1bs=s1bs.process("discard_result", arguments={"data": s1bs})
#job = s1bs.execute_batch(out_format="NetCDF")
s1bs.download("out",format="NetCDF")

If a batch job was used, we can get a description of the job and check its status.

job = connection.job('ea8a0263-9697-4357-a053-db9b6b16e270')
job.describe_job()

{'created': '2021-03-10T08:38:23Z',
 'id': 'ea8a0263-9697-4357-a053-db9b6b16e270',
 'process': {'process_graph': {'ardnormalizedradarbackscatter1': {'arguments': {'data': {'from_node': 'loadcollection1'},
     'elevation_model': None,
     'ellipsoid_incidence_angle': False,
     'noise_removal': True},
    'process_id': 'ard_normalized_radar_backscatter'},
   'loadcollection1': {'arguments': {'bands': ['VV', 'VH'],
     'id': 'SENTINEL1_GRD',
     'spatial_extent': {'crs': 'EPSG:4326',
      'east': 11.469382981297233,
      'north': 46.51363371498709,
      'south': 46.46635946629294,
      'west': 11.218585709568718},
     'temporal_extent': ['2021-01-01', '2021-01-10']},
    'process_id': 'load_collection'},
   'saveresult1': {'arguments': {'data': {'from_node': 'ardnormalizedradarbackscatter1'},
     'format': 'NetCDF',
     'options': {}},
    'process_id': 'save_result',
    'result': True}}},
 'status': 'error',
 'updated': '2021-03-10T08:41:25Z'}

log = job.logs()[0]
log.message

results = job.get_results()
results.get_metadata()

{'assets': {'out': {'card4l:nodata': 0.0,
   'eo:bands': [{'center_wavelength': None, 'name': 'VV'},
    {'center_wavelength': None, 'name': 'VH'},
    {'center_wavelength': None, 'name': 'mask'},
    {'center_wavelength': None, 'name': 'local_incidence_angle'}],
   'href': 'https://openeo-dev.vito.be/openeo/1.0/jobs/65a2e668-b07c-4079-a97f-47012f7ad62c/results/out',
   'type': 'image/tiff; application=geotiff'},
  's1_rtc_02F8D2_N46E011_2021_01_02_MULTIBAND.tif': {'href': 'https://openeo-dev.vito.be/openeo/1.0/jobs/65a2e668-b07c-4079-a97f-47012f7ad62c/results/s1_rtc_02F8D2_N46E011_2021_01_02_MULTIBAND.tif',
   'type': 'image/tiff; application=geotiff'},
  's1_rtc_02F8D2_N46E011_2021_01_02_metadata.json': {'href': 'https://openeo-dev.vito.be/openeo/1.0/jobs/65a2e668-b07c-4079-a97f-47012f7ad62c/results/s1_rtc_02F8D2_N46E011_2021_01_02_metadata.json',
   'type': 'application/json'},
  's1_rtc_02F97C_N46E011_2021_01_03_MULTIBAND.tif': {'href': 'https://openeo-dev.vito.be/openeo/1.0/jobs/65a2e668-b07c-4079-a97f-47012f7ad62c/results/s1_rtc_02F97C_N46E011_2021_01_03_MULTIBAND.tif',
   'type': 'image/tiff; application=geotiff'},
  's1_rtc_02F97C_N46E011_2021_01_03_metadata.json': {'href': 'https://openeo-dev.vito.be/openeo/1.0/jobs/65a2e668-b07c-4079-a97f-47012f7ad62c/results/s1_rtc_02F97C_N46E011_2021_01_03_metadata.json',
   'type': 'application/json'},
  's1_rtc_02FB1E_N46E011_2021_01_07_MULTIBAND.tif': {'href': 'https://openeo-dev.vito.be/openeo/1.0/jobs/65a2e668-b07c-4079-a97f-47012f7ad62c/results/s1_rtc_02FB1E_N46E011_2021_01_07_MULTIBAND.tif',
   'type': 'image/tiff; application=geotiff'},
  's1_rtc_02FB1E_N46E011_2021_01_07_metadata.json': {'href': 'https://openeo-dev.vito.be/openeo/1.0/jobs/65a2e668-b07c-4079-a97f-47012f7ad62c/results/s1_rtc_02FB1E_N46E011_2021_01_07_metadata.json',
   'type': 'application/json'},
  's1_rtc_0435A9_N46E011_2021_01_01_MULTIBAND.tif': {'href': 'https://openeo-dev.vito.be/openeo/1.0/jobs/65a2e668-b07c-4079-a97f-47012f7ad62c/results/s1_rtc_0435A9_N46E011_2021_01_01_MULTIBAND.tif',
   'type': 'image/tiff; application=geotiff'},
  's1_rtc_0435A9_N46E011_2021_01_01_metadata.json': {'href': 'https://openeo-dev.vito.be/openeo/1.0/jobs/65a2e668-b07c-4079-a97f-47012f7ad62c/results/s1_rtc_0435A9_N46E011_2021_01_01_metadata.json',
   'type': 'application/json'},
  's1_rtc_04393F_N46E011_2021_01_08_MULTIBAND.tif': {'href': 'https://openeo-dev.vito.be/openeo/1.0/jobs/65a2e668-b07c-4079-a97f-47012f7ad62c/results/s1_rtc_04393F_N46E011_2021_01_08_MULTIBAND.tif',
   'type': 'image/tiff; application=geotiff'},
  's1_rtc_04393F_N46E011_2021_01_08_metadata.json': {'href': 'https://openeo-dev.vito.be/openeo/1.0/jobs/65a2e668-b07c-4079-a97f-47012f7ad62c/results/s1_rtc_04393F_N46E011_2021_01_08_metadata.json',
   'type': 'application/json'},
  's1_rtc_043A0E_N46E011_2021_01_09_MULTIBAND.tif': {'href': 'https://openeo-dev.vito.be/openeo/1.0/jobs/65a2e668-b07c-4079-a97f-47012f7ad62c/results/s1_rtc_043A0E_N46E011_2021_01_09_MULTIBAND.tif',
   'type': 'image/tiff; application=geotiff'},
  's1_rtc_043A0E_N46E011_2021_01_09_metadata.json': {'href': 'https://openeo-dev.vito.be/openeo/1.0/jobs/65a2e668-b07c-4079-a97f-47012f7ad62c/results/s1_rtc_043A0E_N46E011_2021_01_09_metadata.json',
   'type': 'application/json'}},
 'bbox': [11.218585709568718,
  46.46635946629294,
  11.469382981297233,
  46.51363371498709],
 'geometry': {'coordinates': [[[11.218585709568718, 46.46635946629294],
    [11.218585709568718, 46.51363371498709],
    [11.469382981297233, 46.51363371498709],
    [11.469382981297233, 46.46635946629294],
    [11.218585709568718, 46.46635946629294]]],
  'type': 'Polygon'},
 'id': '65a2e668-b07c-4079-a97f-47012f7ad62c',
 'links': [{'href': 'https://openeo-dev.vito.be/openeo/1.0/jobs/65a2e668-b07c-4079-a97f-47012f7ad62c/results',
   'rel': 'self',
   'type': 'application/json'},
  {'href': 'http://ceos.org/ard/files/PFS/SR/v5.0/CARD4L_Product_Family_Specification_Surface_Reflectance-v5.0.pdf',
   'rel': 'card4l-document',
   'type': 'application/pdf'}],
 'properties': {'created': '2021-03-10T07:30:07Z',
  'datetime': None,
  'end_datetime': '2021-01-10T00:00:00Z',
  'instruments': ['c-sar'],
  'processing:lineage': {'process_graph': {'ardnormalizedradarbackscatter1': {'arguments': {'data': {'from_node': 'loadcollection1'},
      'elevation_model': None,
      'ellipsoid_incidence_angle': False,
      'noise_removal': True},
     'process_id': 'ard_normalized_radar_backscatter'},
    'loadcollection1': {'arguments': {'bands': ['VV', 'VH'],
      'id': 'SENTINEL1_GRD',
      'spatial_extent': {'crs': 'EPSG:4326',
       'east': 11.469382981297233,
       'north': 46.51363371498709,
       'south': 46.46635946629294,
       'west': 11.218585709568718},
      'temporal_extent': ['2021-01-01', '2021-01-10']},
     'process_id': 'load_collection'},
    'saveresult1': {'arguments': {'data': {'from_node': 'ardnormalizedradarbackscatter1'},
      'format': 'GTiff',
      'options': {}},
     'process_id': 'save_result',
     'result': True}}},
  'proj:epsg': 3785,
  'start_datetime': '2021-01-01T00:00:00Z',
  'updated': '2021-03-10T07:36:34Z'},
 'stac_extensions': ['processing', 'card4l-eo', 'eo', 'projection'],
 'stac_version': '0.9.0',
 'type': 'Feature'}

Download results

CARD4L results contain STAC metadata and the requested image.

Here we simply download everything, for inspection.

results.download_files()

[PosixPath('/data/users/Public/driesj/openeo/SRR1_notebooks/out'),
 PosixPath('/data/users/Public/driesj/openeo/SRR1_notebooks/s1_rtc_02F8D2_N46E011_2021_01_02_MULTIBAND.tif'),
 PosixPath('/data/users/Public/driesj/openeo/SRR1_notebooks/s1_rtc_02F8D2_N46E011_2021_01_02_metadata.json'),
 PosixPath('/data/users/Public/driesj/openeo/SRR1_notebooks/s1_rtc_02F97C_N46E011_2021_01_03_MULTIBAND.tif'),
 PosixPath('/data/users/Public/driesj/openeo/SRR1_notebooks/s1_rtc_02F97C_N46E011_2021_01_03_metadata.json'),
 PosixPath('/data/users/Public/driesj/openeo/SRR1_notebooks/s1_rtc_02FB1E_N46E011_2021_01_07_MULTIBAND.tif'),
 PosixPath('/data/users/Public/driesj/openeo/SRR1_notebooks/s1_rtc_02FB1E_N46E011_2021_01_07_metadata.json'),
 PosixPath('/data/users/Public/driesj/openeo/SRR1_notebooks/s1_rtc_0435A9_N46E011_2021_01_01_MULTIBAND.tif'),
 PosixPath('/data/users/Public/driesj/openeo/SRR1_notebooks/s1_rtc_0435A9_N46E011_2021_01_01_metadata.json'),
 PosixPath('/data/users/Public/driesj/openeo/SRR1_notebooks/s1_rtc_04393F_N46E011_2021_01_08_MULTIBAND.tif'),
 PosixPath('/data/users/Public/driesj/openeo/SRR1_notebooks/s1_rtc_04393F_N46E011_2021_01_08_metadata.json'),
 PosixPath('/data/users/Public/driesj/openeo/SRR1_notebooks/s1_rtc_043A0E_N46E011_2021_01_09_MULTIBAND.tif'),
 PosixPath('/data/users/Public/driesj/openeo/SRR1_notebooks/s1_rtc_043A0E_N46E011_2021_01_09_metadata.json')]

import hvplot.xarray
sar = xr.open_dataset('out',engine="h5netcdf")
sar

<xarray.Dataset>
Dimensions:                (t: 6, x: 1684, y: 577)
Coordinates:
  * t                      (t) datetime64[ns] 2021-01-01 ... 2021-01-09
  * x                      (x) float64 6.717e+05 6.717e+05 ... 6.885e+05
  * y                      (y) float64 5.149e+06 5.149e+06 ... 5.154e+06
Data variables:
    VH                     (t, y, x) float32 ...
    VV                     (t, y, x) float32 ...
    local_incidence_angle  (t, y, x) float32 ...
    mask                   (t, y, x) float32 ...
Attributes:
    crs:      +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
    nodata:   0.0

import numpy as np
(10* np.log10(sar.VH)).hvplot.hist( bin_range=(-30, 0))

sar.t

<xarray.DataArray 't' (t: 6)>
array(['2021-01-01T00:00:00.000000000', '2021-01-02T00:00:00.000000000',
       '2021-01-03T00:00:00.000000000', '2021-01-07T00:00:00.000000000',
       '2021-01-08T00:00:00.000000000', '2021-01-09T00:00:00.000000000'],
      dtype='datetime64[ns]')
Coordinates:
  * t        (t) datetime64[ns] 2021-01-01 2021-01-02 ... 2021-01-08 2021-01-09

import holoviews as hv
(10* np.log10(sar.VH)).sel(t="2021-01-02").hvplot(cmap='RdYlBu',dynamic=False,width=1300,height=1300).redim(VH=hv.Dimension('VH', range=(-25, -5)))

Timeseries comparison

Next to requesting CARD4L backscatter data on SentinelHub, we can also generate Sigma0 backscatter on CreoDIAS.

backend_url = "https://openeo.creo.vito.be"
con = openeo.connect(backend_url).authenticate_basic()

W, S = 510000, 5680000
bbox = {
    "west": W, "east": W + 256 * 10,
    "south": S, "north": S + 256 * 10,
    "crs": 32631
}

dates = ("2020-05-06T00:00:00", "2020-09-30T00:00:00")

def backscatter(connection,noise=True):
    return (connection.load_collection("SENTINEL1_GRD")
    .filter_bbox(**bbox)
    .filter_temporal(dates)
    .filter_bands(["VH", "VV"])
    .sar_backscatter(coefficient="sigma0-ellipsoid", noise_removal=noise, options={"tile_size": 256})
    .apply(lambda x: 10 * x.log(base=10)))

%%time
backscatter(con,noise=True).download("creo-series.nc", format="NetCDF")
backscatter(con,noise=False).download("creo-series-noise.nc", format="NetCDF")

CPU times: user 1.96 s, sys: 459 ms, total: 2.41 s
Wall time: 16min 30s

%%time
backscatter(connection).download("shub-series.nc", format="NetCDF")

CPU times: user 560 ms, sys: 274 ms, total: 834 ms
Wall time: 55.8 s

%%time
asc = (connection.load_collection("S1_GRD_SIGMA0_ASCENDING").filter_bbox(**bbox)
    .filter_temporal(dates)
    .filter_bands(["VH", "VV"]))
desc = (connection.load_collection("S1_GRD_SIGMA0_DESCENDING").filter_bbox(**bbox)
    .filter_temporal(dates)
    .filter_bands(["VH", "VV"]))
desc.merge_cubes(asc,overlap_resolver="max").apply(lambda x: 10 * x.log(base=10)).download("snap-series.nc", format="NetCDF") 

CPU times: user 750 ms, sys: 313 ms, total: 1.06 s
Wall time: 48.8 s

creo_ts = xr.open_dataset("creo-series.nc",engine="h5netcdf")
creo_ts

<xarray.Dataset>
Dimensions:  (t: 98, x: 256, y: 256)
Coordinates:
  * t        (t) datetime64[ns] 2020-05-07 2020-05-08 ... 2020-09-28 2020-09-29
  * x        (x) float64 5.1e+05 5.1e+05 5.1e+05 ... 5.125e+05 5.126e+05
  * y        (y) float64 5.68e+06 5.68e+06 5.68e+06 ... 5.683e+06 5.683e+06
Data variables:
    VH       (t, y, x) float32 ...
    VV       (t, y, x) float32 ...
Attributes:
    crs:      +proj=utm +zone=31 +datum=WGS84 +units=m +no_defs 
    nodata:   nan

shub_ts = xr.open_dataset("shub-series.nc",engine="h5netcdf")
shub_ts

<xarray.Dataset>
Dimensions:  (t: 98, x: 256, y: 256)
Coordinates:
  * t        (t) datetime64[ns] 2020-05-07 2020-05-08 ... 2020-09-28 2020-09-29
  * x        (x) float64 5.1e+05 5.1e+05 5.1e+05 ... 5.125e+05 5.126e+05
  * y        (y) float64 5.68e+06 5.68e+06 5.68e+06 ... 5.683e+06 5.683e+06
Data variables:
    VH       (t, y, x) float32 ...
    VV       (t, y, x) float32 ...
Attributes:
    crs:      +proj=utm +zone=31 +datum=WGS84 +units=m +no_defs 
    nodata:   nan

snap_ts = xr.open_dataset("snap-series.nc",engine="h5netcdf")
snap_ts

<xarray.Dataset>
Dimensions:  (t: 98, x: 256, y: 256)
Coordinates:
  * t        (t) datetime64[ns] 2020-05-07 2020-05-08 ... 2020-09-28 2020-09-29
  * x        (x) float64 5.1e+05 5.1e+05 5.1e+05 ... 5.125e+05 5.126e+05
  * y        (y) float64 5.68e+06 5.68e+06 5.68e+06 ... 5.683e+06 5.683e+06
Data variables:
    VH       (t, y, x) float32 ...
    VV       (t, y, x) float32 ...
Attributes:
    crs:      +proj=utm +zone=31 +datum=WGS84 +units=m +no_defs 
    nodata:   nan

Compute the mean value, for each timestamp. OpenEO can also do this for you.

xr.merge([creo_ts, shub_ts.rename({"VV": "VV_SHUB", "VH":"VH_SHUB"}), snap_ts.rename({"VV": "VV_SNAP", "VH":"VH_SNAP"})]).mean(dim=['x','y']).hvplot(width=1200,height=600)